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(57) Abstract 

An automated method and apparatus determines characteristics of nerve fibers. The apparatus includes an imaging unit 
that creates a digital image data representative of an image of a cross-section of a nerve fiber to be tested; a control unit, includ- 
ing image processing capability, is coupled to the imaging unit and analyzes the digital image data created by the imaging unit to 
identify particular characteristics of the nerve fibers including axon area, the myelin sheath area, the fiber area, the myelin sheath 
thickness, the fiber perimeter and the fiber eccentricity. 
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AUTOMATED METHOD AND APPARATUS FOR 
DETERMINING CHARACTERISTICS OF NERVE FIBERS 



Background of the Invention 
The present invention is directed to an automated 
method and apparatus for evaluating fibers in a 
neurological system. In particular, the invention is 
5 directed to an automated method and apparatus that 
employs image analysis techniques to identify and count 
the number of myelinated nerve fibers in a nerve cross- 
section. 

The problem of studying nerve fiber cross sections 
10 quantitatively arises frequently in the' field of 
neuroscience. Parameters such as axon diameter, myelin 
sheath thickness, eccentricity and nerve fiber density 
are useful in many different research areas. 
Measurements of axon diameters and myelin sheath 
15 thicknesses have proven to be of great importance in 
developing mathematical models to describe the 
conduction properties of myelinated nerve fibers. The 
conduction velocity of nerve impulses in a nerve fiber 
can be estimated from the axon diameter, the myelin 
20 sheath thickness, and the internode distances. 

By performing such measurements on a statistically 
significant number of nerve fibers in a nerve bundle, 
and analyzing the distributions of fiber size and 
shape, neurologists can also draw conclusions about the 
• 25 anatomy and physiology of the nervous system. This 

information is significant in assessing the damage to a 
nerve due to natural or experimental causes. 
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The diameters of the axons and the thicknesses of 
the myelin sheaths have Dreviousiy been measured using 
manual techniques and the results entered into a 
computer for processing and storage. Manual 
5 measurement, however, is a very time-consuming process 
because such analyses involves a large n umb er of nerve 
fibers, usually in the range of 1,000 to 3,000 for each 
transverse nerve section. Human errors are unavoidable 
in such a laborious task. Moreover, manual observation 
10 of objects seen under a microscope often calls for 
subjective judgments and decisions which causes 
inconsistent results when several microscope operators 
are employed on a single project. 

Attempts have been made in the past to automate 
15 the process described above. For example, Belson, 
Dudley and Ledley reported in 1968 the use of computer 
to automatically measure the size of neurons. The 
system employed an IBM 7 094 computer equipped with a 
FIDAC scanning instrument. See, "Automatic Computer 
20 Measurements of Neurons" by M. Belson et al., Pattern 
Recognition, 1:119-128, 1986"; and "High-speed 
Automatic Analysis of Biomedical Pictures 11 , by R.S. 
Ledley, Science, 146:216-223, 1964. The system was 
designed to measure the size of the nucleus and the 
25 perikaryon of the neuron, from which information the 
cell volume is derived. 

In 1975, Dunn, O'Leary and Kumley described a 
system developed to count nerve fibers. See, 
"Quantitative Analysis of Micrographs by Computer 
30 Graphics", R.F. Dunn et al., Journal of Microscopy, 
105:205-213, 1975; and "On-line Computerized Entry and 
Display of Nerve Fiber Cross Sections Using Single or 
Segmented Histological Records", Computer and 
Biomedical Research, 9:229-237, 1976. Each fiber in 
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this system was characterized by eight points input 
manually by the operator through a graphics tablet. 
Four of those points were used to indicate the outer 
boundary of the myelin sheath and the remaining four 
5 were used for the boundary of the axon. The computer 
system then searched for the best set of two concentric 
circles to represent that particular fiber. The 
process , however, involved mostly manual tasks and the 
computer was not responsible for the detection of nerve 
10 fibers in an image or the shape analysis of the nerve 
fibers. Instead , the computer merely searched for the 
set of concentric circles which best matched the eight 
points * input by the operator. Thus a mistake made by 
the operator would result in an incorrect output from 
15 the computer. 

A similar method was used by Geckle and Jenkins to 
automatically trace the boundary of the axon. See, 
"Computer-aided Measurements of Transverse Axon 
Sections for Morphology Studies", by W.J. Geckle and 
20 R.E. Jenkins, Computer and Biomedical Research, 16:287- 
299, 1983. The user, however, has to provide an 
arbitrary starting point somewhere inside the axon via 
a digitizing tablet. In other words, the detection of 
the nerve fibers is manual, while the shape of nerve 
25 fibers is traced automatically. The algorithm used is 
a variation to the blob counting method described in 
Computer Vision , by D.H. Ballard and CM. Brown, 
Prentice-Hall, Englewood Cliffs, New Jersey.. 

A more sophisticated system has been developed by 
30 Ellis, Rosen and Cavanagh. See, "Automated Measurement 
of Peripheral Nerve Fibers in Transverse Section", by 
T.J. Ellis et al., Journal of Biomedical Engineering, 
2:272-280, 1980? and "A Programmable Flying-spot 
Microscope and Picture Pre-processor " , by M.J. Eccles 
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et al., Journal of Microscopy, 106:33-42, 1976. This 
system attempts to locate the nerve fibers by the blob 
counting algorithm, then lit examines each located 
object to determine whether or not it is . a conjugated 
5 set of nerve fibers. If it is, the object is 
segmented. On the other hand, if the object is 
determined to be a single fiber or an inter-fiber hole, 
a circularity test decides whether it should be 
classified as a fiber, based on the assumption that 
10 axons in general are circular and that inter-fiber 
holes usually contain irregular protrusions and can be 
rejected by testing for circularity. 

The different schemes described above were 
designed to deal with one very special category of 
15 histological slides - those containing circular 
objects. A certain degree of automation is achieved in 
each method, but all of them require manual inputs of 
some sort which is undesirable for the reasons 
previously stated. Full automation of the process, 
20 however, presents a very difficult set of problems that 
have not been addressed by the conventional methods 
described above. 

In general, the procedure of analyzing the 
histological slide images can be broken down into two 
25 parts. First, the objects have to be identified. This 
is difficult to automate because the computer has to be 
programmed to understand an image which is no more than 
a very large set of pixels with different intensity 
values, requiring immense computational effort and 
30 capacity. Second, morphological data are obtained from 
the identified objects, a step which is very dependent 
on the success of the identification step. 

One of the similarities in most of the previous 
attempts is to start with a manual thresholding, which 
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converts a gray level image into a one-bit image - an 
image which consists only of boundary (myelin sheath) 
and non-boundary (axon or intercellular space) pixels. 
The advantage in using manual thresholding is that 
5 analyzing one-bit images is much easier than analyzing 
gray scale images since the axon boundaries in 1-bit 
images are clearly defined. However, manual 
thresholding has its own shortcomings. Namely, the 
results of manual thresholding are undesirably 
10 inconsistent and not reproducible in general. 

The detection of the objects in all the previous 
works described requires successful thresholdings. The 
objects can easily be detected by using blob counting 
if they are nicely separated from each other. This is 
15 usually not the case in practice, however, as the 
images of the nerve fibers often appear to be clumped 
together due to the high packing density of the nerve 
fibers and the limited resolving power of the light 
microscope. A segmentation algorithm must therefore be 
20 employed as each cluster will be identified as one 
object using blob counting. 

The above described situation can be improved by 
using electron micrographs instead of light 
micrographs. See, "Morphometric Studies on the Studies 
25 on the Relationship Between Myelin Sheath and Axon Area 
in the Sciatic Nerve of Adult Rats", by J Ashami et 
al., Anatomischer Anzeiger, 140:52-61, 1976. A higher 
magnification, however, requires a larger number of 
slides to cover the same region and therefore decreases 
30 the efficiency of the system. The location of the 
objects also turns out to be quite complicated and 
time-consuming . 

Of the four groups described above, only Ellis et 
al. have attempted to solve this problem. Their 
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algorithm segments an object based on the local 
curvature aiong the boundary. The local curvature is 
first smoothed by a convolution procedure as described 
in "Analysis of Digitized Boundaries of Planar 
5 Objects:, by M.J. Eccles et al., Pattern Recognition, 
9:31-41, 1977. Turning points of local curvatures are 
identified as nodes. Nodes are matched to specify 
lines along which the object is segmented. The 
matching is dependent on the orientation, and the 
10 absolute distance between, the nodes. 

There are several disadvantages associated with 
the above-described procedure. For example, it is not 
clear how many smoothing convolutions are needed. 
Further, the segmentation algorithm assumes the 
15 boundaries of the objects to be continuous. The 
segmentation algorithm, however, is employed after a 
single level manual thresholding that has been found to 
produce a significant amount of discontinuities in the 
object boundaries. Measurement of identified objects 
20 is also the last part of the automatic analysis which 
involves boundary tracing that traditionally employs 
the blob analysis. Obviously, if segmentation were 
unsatisfactory, the blob counting that followed would 
be unsuccessful. 
25 It is therefore an object of the present invention 

to overcome the deficiencies of the prior art methods 
and systems discussed above by providing a system that 
is capable of digitizing images of histological slides 
from the cross-section of cochlear nerves, counting the 
30 number of myelinated nerve fibers in each slide, 
obtaining morphological statistics about the nerve 
fibers (including axon area, myelin sheath thickness, 
fiber eccentricity, etc.) and assessing the damage of 
cochlear nerves due to natural or experimental causes , 
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including damages induced by electrical stimulation as 
in cochlear implants . 

A further object of the invention is to provide a 
system that can employ a relatively low cost, general 
5 purpose computer and use the lowest possible 
magnification in order to include as many nerve fibers 
in one frame as possible, so that the number of frames 
needed to be processed will be reduced to a minimum. 

A still further object of the invention is to 
10 provide a system that can accomplish the above-noted 
function at a reasonable speed and accuracy, and with a 
reduction in the amount of human inputs . 

Summary of the Invention 
A method and apparatus is disclosed for 

15 determining characteristics of nerve fibers. The 
apparatus includes an imaging unit that creates a 
digital representation of an image of a crbss-section 
of a nerve fiber to be tested, and a control unit. The 
control unit includes image processing capability, and 

20 is coupled to the imaging unit. The control unit 
analyzes the digital image data created by the imaging 
unit to identify particular characteristics of the 
nerve fiber. 

In operation, digital image data representative of 
25 an image of a nerve fiber cross-section are generated 
by the imaging unit, for example an electronic camera, 
and are stored for processing by the control unit. The 
control unit accesses the stored digital image data to 
locate local maxima of the digital image data, performs 
30 a variable thresholding operation, refines the local 
maxima, defines a plurality of blob image regions using 
blob coloring, and performs an elimination operation to 
identify which of the blob image regions represent 
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nerve fibers. The control unit further calculates at 
least one of the axon area, the myelin sheath area, the 
fiber area, the myelin sheath thickness , the fiber 
perimeter and the fiber eccentricity of the blob image 
5 regions that represent nerve fibers . A unique analysis 
based on eccentricity characteristics is used to 
eliminate blob regions that are not axons . 

Brief Description of the Drawings 
With the above as background, reference should now 
10 be made to the following detailed description of the 
invention with reference to the accompanying drawings, 
in which: 

Fig. 1.1 illustrates a cross-section of a normal 
cochlear nerve from a guinea pig; 
15 Fig. 1.2 illustrates a cross-section of a cochlear 

nerve from a guinea pig after electrical stimulation; 

Fig. 2.1 illustrates a block diagram of a system 
in accordance with the present invention; 

Fig. 2.2 illustrates a block diagram of a FG-100- 
20 AT-1024 image processor board employed in the system 
illustrated in Fig. 2.1; 

Figs. 3. 1(a) -(e) illustrate histogram equalization 
wherein Fig. 3.1(a) represents an original image, Fig. 
3.1(b) its histogram, Fig. 3.1(c) a processed picture, 
25 Fig. 3.1(d) its histogram and Fig. 3.1(e) a mapping 
function ? 

Fig. 4.1 shows the 4-connected and 8-connected 
neighbors of a point P; 

Fig. 4.2 shows an example of dilation; 
30 Fig. 4.3 shows an example of erosion; 

Fig 5.1 illustrates the selection of a region in a 
MANUAL operation mode; 
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Fig. 5.2 illustrates the selection of a fiber in 
an AIDED mode of operation? 

Figs. 5.3(a)-(f) illustrate a process of 
disconnecting a fiber from its neighbors and smoothing 
5 a fiber boundary; 

Fig. 5.4 illustrates obtained local maxima? 
Fig. 5.5 illustrates a 16 x 12 image region? 
Fig. 5.6 illustrates a thresholded image? 
Figs. 5. 7 (a) -(b) represent a model of connection 
10 between fibers and inter-fiber holes? 

Fig. 5.8 shows an image of a nerve section after 
blob coloring? 

Fig. 5.9 shows the completed result of a nerve 
section counted by the AUTO mode of operation; 
15 Fig. 6.1 is a table which illustrates the results 

obtained from the system illustrated in Fig. 2.1 in the 
MANUAL, AIDED and AUTO modes of operation? 

Fig. 6.2 is a table which illustrates the analyses 
of six histological sections using the AUTO mode of 
20 operation — three from normal nerve and three from 
electrically stimulated nerves? and 

Fig. 6.3(a)-(f) illustrate six parameters 
generated by the AUTO mode of operation for the image 
illustrated in Fig. 1.2. 

25 

Detailed Description of the Invention 
Additional background information will be 
presented prior to discussing the specific structural 
and functional details of the invention, in order to 
30 provide additional basis for understanding the 
operation of the invention. The background material 
relates to the desirability of producing a system 
capable of identifying degeneration of the cochlear 
nerve. It should be understood, however, that the 
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present invention is not limited to this particular 
application, but is applicable broadly to many 
applications, including, fJr example, peripheral and 
central nerve cross-sections, cytology automation, red 
and white blood cell counts, muscle fiber analysis, 
cell culture colony counting and analysis of an image 
consisting of closely-packed circular objects. 

1.0 Background 

As stated above, the problem of studying nerve 
fiber cross sections quantitatively arises frequently 
in the field of neuroscience. For example, 
neurobiologists have long been interested in the size 
distribution and the regional differences of the 
cochlear nerve fibers . The neurobiologists associated 
the degeneration of the cochlear nerve with 
morphological changes in the nerve fibers as well as 
the changes in the size of the nerve fiber population. 

One of the reasons for studying the cochlear nerve 
in particular stems from the recent rapid developments 
in cochlear implants. The possibility of the 
perception of sound became a reality for profoundly 
deaf individuals with the approval by the FDA of the 
first cochlear implant device in November 1984. Even 
when hearing is lost, some spiral ganglion cells may 
remain which can be electrically stimulated to produce 
the sensation of hearing. 

It has been found, however, that chronic 
electrical stimulation may result in substantial tissue 
damage. See, for example, "Mechanisms of Electrically 
Induced Damage After Cochlear Implantation", Annals of 
Otology, Rhinology and Laryngology, 95:185-189, 1986. 
As successful implants require healthy spiral ganglion 
cells and cochlear nerve fibers, a measure of merit for 
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cochlear prostheses wouid be an assessment of the 
degree of cochlear nerve survival after simulation. 
One way to do this is would be to count the number of 
surviving cochlear nerve fibers from the cross- 
5 sectional view of the cochlear nerve. 

The difference between normal cochlear nerve 
fibers and damaged cochlear nerve fibers will be 
illustrated with reference to two sets of histological 
slides analyzed. One set contains transverse cross- 
10 sections of cochlear nerve from normal guinea pigs. 
The other set is from guinea pigs which have undergone 
a variable amount of cochlear electrical stimulation. 
Two-micron thick sections are obtained every 100 
microns along the cochlear nerve in such a way that the 
15 cochlear nerve fibers appear as circular as possible. 

A cross-section of a cochlear nerve from a guinea 
pig is illustrated in Figure 1.1. The myelin sheaths 
of the nerve fibers appear dark when stained with 
toluidine blue. The nerve fiber diameters range from 
2 0 2-4 jim with an average of about 3jim. The nerve fibers 
are mostly circular and are innervated by blood vessels 
which appear as circular but significantly large 
structure. 

A similar cross-section from a cochlear nerve of a 
25 guinea pig after electrical simulation is shown in 
Figure 1.2. Abnormal crescent-shaped fibers can be 
seen. The presence of these irregular fibers is may be 
due to the activities of microglia - an assortment of 
phagocytic cells that are normally at rest. Under 
30 pathological conditions they change to reactive 
microglia and form processes which pinch into the 
neurons . 
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2 . 0 System Structure 

Referring now to Fig. 2.1, a block diagram of a 
system capable of accurately and automatically 
identifying healthy nerve fibers from damaged nerve 
5 fibers quickly and economically is disclosed. The 
system includes a microscope 10, a video camera 12 , a 
general purpose computer 14 equipped with a special 
purpose image processing board, and an output device 16 
such as graphics monitor or printer. Each of the 
10 different components will be discussed in greater 
detail below. 

2.1 Microscope 

The analysis of cochlear nerve fibers is a tedious 
task because of the large number of miniature fibers 
15 present. For example, the average number of cochlear 
nerve fibers in the guinea pig is as many as 24,000, 
most of which range from 2-4um in diameter. Electron 
microscopes are often used in this kind of analysis in 
order to gain the required resolution. For example, it 
20 is necessary to be able to resolve to about one eighth 
of the diameter of a fiber. That means the image 
should be magnified to such a degree that about 8 
pixels span the diameter of each fiber. A decrease in 
resolution may blur the image and present problem to 
25 any computer algorithm that attempts to extract 
features from the image. An increase in resolution 
which can only be achieved by using the electron 
microscope, however, results in a much higher 
processing time and reduces the number of nerve fibers 
30 that can be processed per frame. 

Since the objects being considered are very small, 
it is important to obtain the necessary resolution. 
The resolution acquirable by a compound microscope is 
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(2.1) 

where X is the wavelength of the light used and N 
is the numerical aperture of the objective: 

/V = n sin ti (2.2) 



5 where u is the half angular aperture of the 

objective and n is the refractive index of the medium. 

If it is assumed that the smallest .fiber 
considered is 2nm in diameter, then the resolution 
required will be 

10 8 (2#3) 

= 0.25am 

If white light is utilized, the applicable wavelength 
would be for green (550nm), which is at the middle of 
the visible spectrum and the region of highest acuity. 
The numerical aperture required according to -Equation 
15 (2.1) is 



v _ °- 61 ^ 0.61 x 550nm 



From the above it is therefore obvious that the 
required resolution will not be obtained by using dry 
20 lenses, because N fiir is one and from Equation (2.2) it 
is known that N a . must be smaller than one. Thus, an 

all 

immersion objective used with immersion oil 
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(n Qii =1.515) or liquid with a higher refractive index 
must be used- For example, a Zeiss Axioplan microscope 
with a plane apochromatic immersion objective having a 
magnification of 63X and a numerical aperture of 1.4 
5 can be employed for the microscope 10 illustrated in 
Fig. 2.1. This is a high-NA objective corrected for 
spherical aberration in four wavelengths, chromatic 
aberration for more than four wavelengths and for 
flatness of fieid. 

10 2.2 Camera 

The camera 12 is directly connected to the 
microscope 10. The camera 12 is preferably a Hitachi 
KP-C100 camera which incorporates a CCD (Charge- 
Coupled Device) solid state detector. The CCD operates 
15 by transferring packets of minority charge from one 
potential well to another. The CCD has a two- 
dimensional grid of sensor elements made of silicon 
photodiodes. The electrons produced in each photodiode 
are immediately transferred to an adjoining, 
20 electrically positive potential well. The electron 
packets trapped in the potential wells are shifted 
upwards row by row, such that the top most row falls 
into a structure called the isolation transfer gate. 
Within the isolation transfer gate, the packets are 
25 shifted to the left one by one and amplified to produce 
an output resembling that of a raster scan system. The 
overall operation of CCD equipped cameras is well known 
in the art. The Hitachi camera produces a picture of 
510x492 pixels. 

30 The scan area of the camera is 8 . 8mmx6 . 6mm. Using 

a 63X objective, the area covered is 
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S800 6600 

UTT. 

(2.5) 



63 ' m * — * m 



= 140um x IOSzxtt: 



The total number of cochlear nerve fibers that can fit 
into this space is about five hundred. 

2.3 Image Processor 
5 The output of the camera 12 is digitized and 

stored for future references and manipulations. This 
task is accomplished by an image processor board 
contained within the general purpose computer 14. 
Alternatively, the image processor can take the form of 
10 a separate unit coupled to the general purpose computer 
14. Preferably, an Imaging Technology imaging board, 
Model FG-100-AT-1024, is used as the image processor, 
although other image processing devices may be readily 
employed. The FG-100-AT-1024 is a single board, real- 
15 time image processor designed to be interfaced with IBM 
PC-AT's and compatibles. A . block diagram of the 
processor is shown in Figure 2.2. The image processor 
digitizes the 512x480 pixel images from the camera with 
8-bit precision ( 256 gray levels). The digitized image 
20 is passed through the feedback/input look-up table and 
stored in a frame memory. The frame memory is arranged 
as a 1024x1024 two-dimensional array of 12-bit pixels. 
It is therefore large enough to store 4 images at the 
same time. . While the lower 8 bits of the pixels are 
25 usually used to store images, the upper 4 bits can be 
used to store other information. 

The input/ feedback look-up table image processor 
takes as input a 12-bit pixel and maps it to another 
12-bit pixel. It serves two functions. It can be used 
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for input preprocessing, and it can also be used 
through a feedback loop to manipulate data already 
present in the frame memory. This look-up table can 
be programmed to perform different tasks in real time. 
5 The output circuitry of the image processor 

involves three look-up tables - one for each primary 
color. These look-up tables provide 12-bit to 8-bit 
mappings and can be programmed for pseudo-color 
processing if desired. 

10 2.4 Computer 

The general purpose computer 14 illustrated in 
Fig. 2.1 acts as a central control unit and also 
controls the operation of the image processor board 
described above. Preferably a Compaq Deskpro 386/25 

15 personal computer (PC) is employed for the general 
purpose computer 14. The frame memory of the image 
processor is memory mapped to the PC which performs 
most of the complicated computations that the image 
processor cannot handle. This includes searching and 

20 other operations requiring random access of the frame 
memory. With careful planning, the image processor can 
perform selected functions in parallel with the 
computer. The use of a PC also facilitates the 
transfer of data to other software packages that are 

25 readily available. The user interacts with the 
computer via a keyboard and a mouse. 

2.5 Output Device 

The output device is provided to allow the user to 
view the data and the results obtained. In addition to 
30 the alphanumeric monitor of the PC, a graphics monitor 
and graphics printer are preferably employed. A 
graphics monitor capable of displaying the entire video 
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frame is required. Because of the degradation of video 
images near the edges, the periphery of a video image 
is usually not shown on commercial TV. A monitor which 
can display the entire frame is said to have the 
5 underscan capability. For example, a Sony KV-1250 
color monitor can be used for the graphics display. 

3 . 0 Image Processing Techniques 

A general description and review of various imaoe 
processing techniques that are employed by the present 
10 invention will be presented prior to discussing the 
specific operational details of the invention. 

3.1 Image Definition 

A digital image will be defined for the purposes 
of describing the invention as a two-dimensional 
15 collection of points with intensity i(x,y) at 
coordinates (x,y). Since intensity is never negative: 

0 < i (x,y) < oo (x , y) 

For the 8-bit images with 512x4 80 pixels provided by 
the camera of the preferred embodiment: 

20 0 < x < X - 1 

0 < y < Y - l 

0 < i(x,y) < 255 

where X = 512 and Y = 480. 

The term "resolution" as applied to -digital image 
15 processing refers to both the sampling frequency in the 
spatial domain and the quantization levels of each 
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pixel. The spatial resolution in the pictures produced 
by the camera is 512x480. The resolution can be 
decreased to 256x240 by replacing every four pixels by 
one pixel having an intensity equal to the average of 
5 the original four in order to remove the redundancies 
when the magnification is too high. It is preferred, 
however, to control the magnification rather than the 
resolution so that each image contains the optimum 
number of fibers . 
10 The quantization level in the pictures produced is 

8-bit , or equivalently 256 gray levels. It is also 
possible, for example, to decrease the quantization 
level to 7 -bits by grouping every two consecutive gray 
levels into one. Decreasing either the spatial 
15 resolution or the quantization levels will blur an 
image but will also decrease the amount of space 
required to store a picture. For example a 512x480, 8- 
bit picture requires 240kb to store, while a 256x240, 
7 -bit picture requires only 5 2.5kb. 
20 Operations on an image can be classified into 

several categories : 

Local Operations - Each pixel in the resulting image 
depends only on its immediate neighborhood - usually 
25 the 3x3 neighborhood - in the original image (e.g. 
spatial domain filtering). 

Regional Operations - Each pixel in the resulting image 
depends on the surrounding region in the original image 
(e.g. regional averaging). 

30 Global Operations - Each pixel in the resulting image 
depends on all the pixels in the original image (e.g. 
Fourier Transform). 
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3.2 Histogram Equalization 

While most image processing techniques may fall 
into any of the above categories, pattern recognition 
techniques are usually local or regional, but rarely 
5 global. 

A histogram of a picture is a plot of intensity 
versus the frequency of occurrence. The range of 
intensity of a picture is often small compared to the 
range available on a system. In this case, it is 

10 possible to improve the contrast by using histogram 
equalization which has the effect of spreading out the 
intensity. An intensity mapping function according to 
which each original gray level is mapped to a new one 
is required to provide equalization. 

15 For the original histogram, the function n (i) is 

the number of pixels that have a gray level i. 
Similarly, the. desired histogram will have the function 
n d (i). The probability density functions (PDF's) are: 

N 

where N is the total number of pixels in the picture, 
20 specifically X x Y = 512x480 = 245760 pixels. In the 
ideal case: 

Mi) = j 

where I = 256 is the total number of gray levels of the 
system. 

Letting the cumulative distribution functions 
25 (CDF's) be then: 

*W -£/.(;) 

i 

m HuU) 
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Since it is impractical to equate f\ with f^, 
instead can be equated with g^: 

9 <(i) = i 



A mapping function M is now needed such that i rf = M 



5 Solving the above equations gives : 

/ 

i<« /V/(i.) s — X n *0) (3.1) 

Each pixel in the original picture can now be 
changed to a new value based on Equation (3.1). 
Histogram equalization is a global process since the 
10 new value of each pixel depends on the histogram of the 
original picture. In the discrete case, it is 
irreversible because two values may be mapped into one. 
Figure 3.1 illustrates an example of histogram 
equalization. 

15 The two-dimensional discrete Fourier transform is: 
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and f(x,y) becomes i(x,y) for an image. For examples 
of 2-D DFT of images see Computer Image Processing and 
Recognition . E.L. Hall, Academic Press, New York, 1979; 
Digital Picture Processing . A. Rosenfeld and A.C. KaJc, 
Volume 1, Academic Press, Orlando Florida, 2nd Edition, 
1982; and Signals and Systems . A.V. Oppenheim et al., 
Prentice-Hall, Englewood Cliffs, New Jersey, 1983. 

3 . 3 Filter Domain Filtering 

To filter an image in the frequency domain 
requires three steps. First, the Fourier transform of 
the image must be obtained by directly applying 
Equation (3.2). Because of the separability property 
of the 2-D DFT, the same result can be obtained by 
applying the one-dimensional fast Fourier transform 
twice, once horizontally and once vertically. Second, 
a multiplication in the frequency domain must be 
performed to obtain: 



where H(u,v) is the filter. For lowpass filtering, 



where D(u,v) = ✓ ( u +v ) and D Q is the cut-off 
frequency. 

Similarly for highpass filtering, 



(3.3) 




(3.4) 




1 for D > D 0 



for D < Do 



(3.5) 
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The final step is to perform the inverse discrete 
Fourier transform: , 

x-\ r-\ 

9(*,y) = £ £ v)r»"*+?> (3.6) 

An algorithm that is used to compute the forward DFT 
5 can also be used to compute the inverse DFT. If the 
input to the algorithm were G*(u,v), the output would 
be g*(x,y) . 

3 . 4 Spatial Domain Filtering 

In image processing, filtering in the spatial 

10 domain is more common than in the frequency domain. In 
the spatial domain, filtering does not require complex 
arithmetic, but the cut-off will not be as sharp as in 
the frequency domain. It is well known that 
multiplications in the frequency domain correspond to 

15 convolutions in the spatial domain, which means that 

g{z,y) = /(=, y) * h{x,y) = £ £ f{ij)h{x - i,y - ;) < 3 ■ 7 ) 

where h(x,y) is the inverse DFT of the filer H(u,v) in 
Equation (3.3). It should be noted that h(x,y) is an 
XxY array, which is too large to use as a convolution 
20 mask. Letting h'(x,y) be an nxn convolution mask which 
is an approximation to h(x,y), then h'(x f y) can be 
obtained by minimizing the error quantity: 

The convolutions are often restricted to local 
operations by using masks having n=3: 
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! M-i.-i) A(0,-i) Ai l, — I j 

i 
j 

A(x,y;= At-I.O) A(0,0) k(l.O) 

A(-l.l) rt(0,l) A(l,l) 



In some applications it is desirable to have n > 3. i n 
most instances this can be achieved by successively 
applying the 3x3 convolution mask. For example, 
h'(x,y)*h'(x,y) produces a 5x5 convolution mask. The 
associativity property of convolution guarantees: 



10 



15 



(i(x,yj . k(z, y;j . A(x, y j = i( x , y) . [A(x, y) . A(x, 



However the left hand side is more efficient in terms 
of the number of multiplications and additions 
required. If „ « 3 the right hand side requires 25 
multiplications and 24 additions for any (x,y) even 
when the 5x5 convolution mask is precalculated . The 
left hand side requires only 18 multiplications and 16 
additions . 

From Equation (3.7): 
rfx.y) = /f«-l.»-lAl.l) + /(« tIr -i ) A ( o.l) + /(« + I,,-nA ( -i. i :. 

/(* - 1, y)A(l, 0) + /(_-, y)A(O.O) + /(, + i, y)A( - 1 . 0) + (3.8) 

ft* - 1, y + i)A(i t -i) + y + 1)A(0 , _ : ) + /(J + lt + 1)A( _. _ 1} 

Modifying Equation (3.2) for n - 3 gives: 



A r (u,v,= i £ £ A(«,y ;e -'H«*ni (3-9) 



=-i jr=— | 



A convolution mask for lowpass filtering is 
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0.1 Rl 0.1 
0.1 1 0.2 0.1 
0.1 0.1 0.1 



(3.10) 



10 



which is in essence a weighted averaging. It is good 
for removing noise but its indiscriminate use may blur 
the image . 

Combing Equation (3.9) and Equation (3.10): 



Ht(u.v) = 



0.011 0.011 0.011 
0.011 0.111 0.011 
0.011 0.011 0.011 



which verifies that h ' is a low pass filter. 



(3.11) 



Note 



that ' is real and even because h 1 ' is real and even . 
A convolution mask for highpass filtering is : 



i -0.25 -0.25 -0.25 

i 

h h (x 9 y) = ; -0.25 3.00 -0.25 

! -0.25 -0.25 -0.25 



(3.12) 



which is good for edge enhancements but it will aiso 
emphasize noise. 

Combing Equation (3.9) and Equation (3.12): 



0J51 0.361 0.361 



0.351 0.111 0.361 j 



OJfil 0.361 0.061 



I 



15 which verifies that h h ' is a highpass filter. 



(3.13) 



Note 



that H^' is real and even because h^' is real and even. 
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3 . 5 Thresholding 

Thresholding is a technique for converting a gray 
level picture into a binary one: 



5 where i fa (x,y) is the resulting binary image and T is 
the threshold level. 

Thresholding is appropriate only when the picture 
in question consists of two distinct regions and the 
histogram is bimodal . In other cases, it is difficult 

10 to decide on the threshold level, particularly by 
automatic methods . The threshold level is manually set 
to the optimal level that preserves most of the 
boundaries of the axons. It has been found, however, 
that non-uniform illumination of the nerve section 

15 under study can severely hinder the success of a global 
thresholding. Thus, as will be described in greater 
detail, it is desirable to apply a variable threshold 
level on local properties . 

4 . 0 Pattern Recognition Techniques 

20 Pattern recognition techniques refer to the 

algorithms used to extract features from an image. The 
image should be fully enhanced before these techniques 
are implemented. It is convenient to adopt a set of 
notations for the attributes of an image and its 

25 components for the purposes of discussing the 
invention . 

•Pixels or points are indicated by capital 
letters: e.g. P, Q. 





for i(z t y) < T 



(3.14) 



• The x,y coordinates and the intensity of the 
30 pixel P are given by P_ x , P and P i respectively. 
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•The Euclidean distance between two pixels is: 

d pq = *'<( p x " Q x> 2 + ( *y " Q y )2) 

4 . 1 Connectivity 

Two adjacent pixeis P and Q are said to be 4- 
5 connected if and only if D pQ = 1 . They are 8-connected 
if and only if l<D pQ </2 . See Figure 4.1 in which the 
4's and 8's denote the 4 -connected and 8-connected 
neighbors of P. The terms 4-neighbors and 8-neighbors 
refer to the points that are 4-connected and 8- 
10 connected respectively. 

The distance between 4-neighbors is always equal 
to 1 unit. However, the distance between 8-neighbors 
can be either 1 or /2 units . In some applications , it 
is more convenient to assume that the distance between 
15 two diagonally connected 8-neighbors is also 1. 

4.2 Edge Detection 

In most images the information is contained in the 
edges — locations where there is an abrupt change in 
intensity. For this reason edge detection is 
20 frequently used. In some applications, however, it has 
been found that edge enhancement may have an adverse 
effect. 

The Gradient of an imaoe is defined to be: 



25 



I" 3ifr.vi 1 



(4.1) 



While the gradient is a measure of the change in 
intensity, the Laplacian is a measure of the rate of 
change in intensity. Since the edges in an image are 
located at points where the second derivative of 
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intensity passes through zero, the Laplacian of an 
image shows the edges . 

4 . 3 Segmentation 

Segmentation refers to the techniques used to 
5 divide a picture into regions with each region 
representing a feature of a considerable portion 
thereof . There are many methods to segment a picture 
including: 

Merging - First consider all the pixels to be in a 
10 different region. Then merge adjacent pixels if they 
pass a homogeneity test. Apply this method repeatedly 
until no more merging occurs . 

Splitting - Unlike merging consider all the pixels 
to be in the same region. Then perform a homogeneity 
test on each of the four triangular subsection of the 
region. Split the subsection into another four 
triangles if the test fails. Repeat until no more 
splitting occurs . 



15 



20 



Split-and-Merge - Even better results can be 
obtained if merging is performed following splitting, 

A major drawback of the conventional segmentation 
algorithms is that they require a tremendous amount of 
computer memory. if r i s the number of different 
regions and N the number of pixels in a picture, the 
25 amount of memory required to store this information 

is «|log 2 R bytes. For 512x480 pictures , at least 537kb 

would be required for merging- Splitting and split- 

and-merge may require less memory- The amount of 
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memory that we just discussed is the minimum required 
to hold the information. j 

An efficient implementation of the above 
algorithms requires a much larger amount of memory. A 
5 reasonably efficient data structure would be to store 
the coordinates of the pixel belonging to the same 
region in a linked list. This would require 8 bytes 
per pixel or 1. 875Mb, Also the recursive methods 
discussed may cause stack overflows. As a result, the 
10 use of these methods is confined only to small regions. 

4 . 4 Thinning 

Thinning refers to the techniques used to convert 
an object into its skeleton while preserving its 
connectivity, symmetry and shape. It is often used to 

15 reduce line images to unit line width which can be 
represented easily. It only makes sense to perform 
this operation when the object is reasonably elongated. 

The minimum t such that S~*~ is empty can be a 
measure of the thickness of an object. The 

20 elongatedness of an object is defined to be: 

E = ^2 (4.8) 
where A is the area, of the object. 
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Thinning algorithms should not be used if E is 
small. Most thinning algorithms repeatedly remove 
border points from the object based on some criteria 
depending on the application until no more removal 
5 occurs . 

4.5 Mathematical Morphology 

Mathematical morphology is an approach to image 
processing which is based on the shape of the objects 
processed. The following discussion is based on the 
concepts presented by Haralick et al. described in 
"Image Analysis Using Mathematical Morphology", r.m. 
Haralick et al., IEEE Transactions on .'Pattern Analysis 
and Machine Intelligence, PAMI 9: 532-550 , 1987, but the 
equations have been reformulated based on Boolean 
arithmetic instead of set arithmetic to facilitate the 
conversion to computer programs. The following logical 
operations will be used: /\ represents OR, V represents 
AND, and represents NOT. 

4.5.1 Dilation 

Dilation is an operation that spreads out the 
marked pixels and has the effect of reducing noise and 
filling small holes. Let x(x,y) and b(x,y) be two 
binary images, then 



10 



IS 
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a\x,y)®b(x,y)= V V a(i, j) A b(x - i,y - j) 

ItsO j'=0 

(4.9) 

wnere cr is tne cmation operator. 

Fig. 4.2 shows an example of dilation. 

4.5.2 Erosion 

Erosion is an operation that contacts the marked 
5 pixels and has the effect of thinning the object and 
expanding holes. Let a(x,y) and b(x,y) be two binary 
images , then 

«(*. V) e *(*, y ) s /\ /\ a(: t ; ) V 6f: -r,;- y; (4-10) 

where is the erosion operator. Figure 4.3 shows an 

10 example of erosion. 

The most common dilation and erosion operations 
have as input an image and a structuring element known 
as the dilation or erosion mask. The shape of the 
structuring element known as the dilation or erosion 

15 mask depends on the application. 
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4.5.3 Opening and Closing 

Dilation and erosion are often performed in pairs. 
Calling the structuring element k(x,y) and the image 
i(x,y) then: 



i(:,y)Ok(x,y) = (i( Sf yj S fc(x, y)) e k(z,y) (4 . u ) 

wnere © is the opening operation, and 

where ® is the closing operation. (4.12) 



10 
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4 . 6 Blob Counting 

Blob counting usually refers to the techniques of 
locating marked objects and obtaining information about 
them. Assume that the pixels in the objects all have 
value 1 and the background pixels all have value 0. 
The technique for locating the objects is well known. 
The image is scanned until an object pixel (which is 
the starting pixel for boundary tracing) is 
encountered. The boundary is traced by turning left 
when reaching an object pixel and turning right 
otherwise until one arrives back at the starting pixel. 
This procedure will trace the object in a clockwise 
manner. It is essential to mark the boundary pixels as 
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visited so that this object will be ignored in the next 
scan. | 

A convenient way to keep track of the shape of an 
object during boundary tracing is to use the chain code 
representation, wherein each time an object pixel is 
reached, it is considered that a movement has occurred 
from the previous object pixel; 

• A move right is coded as 0 and stored in the 

chain 

• A move down is coded as 1 and stored in the 

chain 

• A move left is coded as 2 and stored in the 

chain 

• A move up is coded as 3 and stored in the chain 
It has been found that the method described by 

Ballard and Brown in Computer Vision has a slight flaw 
which may give an incorrect result when part of the 
object is 1 pixel thick. For example when tracing a 
vertical line of i pixel thick, it is possible to 
return to the starting pixel from the other side of the 
line. Thus, one must make sure that one returns to the 
starting pixel on the starting side to consider the 
procedure finished. This can be done by having an 
additional constraint — The tracing is completed if 
a return to the starting pixel occurs and the chain 
code is either 0 or 3. 
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The perimeter and area of the object can be easily 
obtained from the chain code. The perimeter is simply 
the number of elements in the chain. In order to 
calculate the area we start by initializing the area to 
0 and keep track of the y-coordinate as we traverse the 
chain : 

• Chain code 0: subtract the value of y from the 
cumulative area 

• Chain code 1: increment the value of y 

• Chain code 2: add the value of y to the 
cumulative area 

• Chain code 3: decrement the value of y 



4 . 7 Correlation 

The similarity of a some function g(x,y) to f (x,y) 
15 can be checked by defining the dissimilarity to be: ■ 

X-l Y-\ 



10 



vu,a) = ZZ!/(*.»)-«r,y)p 



r=0 t=o 



20 



_ ™J .., x-iy, ,_ JK _, (4.13) 
C0n3l4nt «M 



where S(f,gj is a measure of similarity. 

In order to search for an object g(x,y) with a 
prior knowledge of its shape in f(x,y), the most 
reliable method is to conduct an extensive search for 
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the x and y coordinates that would maximize the cross- 
correlation function: i 



/(*, V) o g(z, y) = r £ f(i t j)g(x + «, y + j) ( « • " ) 

•»0 >s=0 



f *( x rj) would just become f(.i,j) for images. When 
matching different objects in the picture, and since 
g(x,y) is not constant as shown in Equation 4.13 f a 
normalized cross-correlation function must be employed: 



/(*■ v) • •(*, y) = tt l ll&f m (ij)9(x + i*w + j) (4.15) 



4.8 Hough Transform 

10 The Hough transform is another way to locate 

features in an object. Its most common use is to 
detect straight lines, but it can be generalized to 
detect arbitrary shapes. The approach is to tally the 
points in the parameter space. For example, in order 

15 to detect straight lines in a noisy image, one could 
mark the corresponding point in the slope intercept 
space (mc-space) for each object point identified in an 
image based on the equation y = mx + c. The maximum 
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poinrs in the mc-space indicate a possibility broken 
straight line. 

The same idea can be used for the detection of 
circles , only now the parameter space would be three- 

5 dimensional with x q/ y f and r from the circle equation 

2 2 2 

(x-x Q ) - (y - y )* = r . This means that the space 

requirement would be tremendous. If it is assumed 

0<x Q <511, 0<y o <479, and 0<r<255, reserving one byte per 

voxel for the tallies, the space requirement would be 

10 60Mb. 

5 - 0 System Operation 

The above-described image processing techniques 
are employed by the system illustrated in Fig. 2.1 to 
count and extract data from the histological slides. 
15 The procedures involved are listed in Table 5.2. For 
illustrative purposes, the slide shown in Figure 1.1 
will be used unless otherwise stated. The system is 
capable of operating in three modes (MANUAL, AIDED and 
AUT0) f as will be described in greater detail below. 

20 

5.1 Statistics Generated 

The first step implemented by the system of Fig. 
2.1 is the derivation of the statistics regarding axon 
area ( AA ) , myelin sheath thickness ( MT ) and 
25 eccentricity (EC). Statistics about total fiber area 
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Procedure 


Time Required 


Initialization " , 


0.94% 


Preprocessing * 


4.92% 


Locating Local Maxima • 


1.56% 


Computaion of A. B, C, and D \ 


0.76% 


Variable Thresholding * 


5.70% 


Refining Local Maxima " 


0.02% 


Blob Coloring 


49.80% 


Selective Erosion 


8.72% 


Setting up the Adjacency Lists 


2.90% 


Symbolic Constraint Propagation 


1.63% 

— ' 


Location of Axon and Myelin Sheath Area 


23.047o i 



Tu6k. 
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(TA), fiber perimeter (FP), and myelin sheath area (MA) 
are also included. 

Each fiber is represented by a set of axon pixels 
and a set of myelin pixels. The fiber perimeter (FP) 
is traced using blob counting. Its value is calculated 
based on 6 -connectivity, and is therefore incremented 
by the square root of two units for each diagonal move 
and by one unit for each non-diagonal move. The other 
values are calculated as follows: 



Totai Area = TA = AA -r MA 
Myeiin Thickness = MT = , — _ , 



(5.1 1 



rp2 

£.ccenirici:v =s EC = ~_ , 

TA - a - 3 ' 



Eccentricity allows one to differentiate between 
circular and non-circular fibers. The minimum 
eccentricity is that of a perfect circle which is 4n = 
12.6. Other shapes have higher eccentricities. a 
perfect square, for example, has an eccentricity of 16. 
It is possible, however, for the calculated 
eccentricity of an object represented on a digital grid 
to be lower than 4n. it would be expected that the 
eccentricity of an object displayed on a digital grid 
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wouid be higher than that of its analog counterpart 
since a circle will no longer be circular on a digital 
grid. The problem is that even though the area is 
generally accepted to be the number of pixels the 
5 object occupies, the perimeter is not well defined. 
Some define the perimeter to be the number of boundary 
pixels. But this definition depends on whether we 
consider 4-connectivity or 8 -connectivity . The 
definition of perimeter as employed by the present 
10 invention, and as described before , tends to give 
values close to 4tt for circular objects. The 
difficulty of defining eccentricity and perimeter is 
even greater for small objects. 

It is important to note that the myelin sheath 
15 . thickness is calculated from the axon area and myelin 
sheath area assuming that both the axon and the fiber 
are perfectly circular. It should also be noted that 
all of the parameters can be calculated when given the 
axon area, myelin sheath thickness and eccentricity. 
20 Histograms of all these parameters are generated by the 
system of the present invention. 

5.2 MANUAL MODE 

The MANUAL mode of operation of the system is the 
least sophisticated mode of operation, but it gives the 
25 user the maximum amount of flexibility. The user uses 



WO 92/17848 



PCI7US92/02563 



3 8 



10 



a mouse to trace out the axon and the myelin sheath 
boundaries and the computer 14 is utilized to obtain 
the desired statistics. More specifically, the user 
selects a region of the slide displayed on the monitor 
and an enlarged view of the selected region is 
displayed as shown in Fig. 5.1. The user places the 
cursor on an axon boundary, activates the mouse and 
proceeds to trace the axon boundary of a nerve fiber. 
Upon returning to the initial starting point, the 
computer calculates the axon area from the traced 
boundary and colors in the axon area, for example with 
yellow, to indicate to the user the completion of this 
step. The user then traces the myelin sheath boundary 
in a similar fashion and the computer calculates the 
myelin sheath area and colors it green. Fig. 5.1 
illustrates the intermediate step where the user has 
already traced out the axon boundary and is in the 
process of completing the myelin sheath boundary. Four 
other fibers have already been counted.. 

It is readily apparent that the MANUAL mode of 
operation is very time consuming. The user must trace 
the axon and myelin boundary for each fiber. The 
computer 14 simply calculates the axon and myelin 
sheath areas by determining the number of pixels within 
25 the boundaries defined by the user. The MANUAL mode of 
operation, however, does provide the ability to fill 



15 



20 
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the region when given a boundary that is traced either 
clockwise or counterclockwise. The points on the 
boundary are stored in a linked list called boundary 
list and the computer 14 perforins the following 
5 procedures to fill the region: 

1. Eliminate repeated points in boundary list. 

2. For every pair of- points in boundary list, 
retain only the point with the smaller y-coordinate. 
If the points have the same y-coordinate , then keep 

) both and create a duplicate pair, for example: 

(xi.y) - (z a ,y) 

becomes 

(xi.y) - (xi,y) — (x 3l y) — (r 2> y) 



3 . Sort the boundary list in ascending order 
according to the value y x 1 + x where 1 = 512 is the 
horizontal window size. 

4. For every pair of points in boundary list that 
have the same y-coordinate , draw a line between them if 
this was not done for the last pair of points. 
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5.4 AIDED MODE 

The AIDED mode of operation is similar to that 
described by Geckle and Jenkins discussed above. It is 
not totally automatic but requires an operator to 
5 identify the fibers. In order to count a fiber, the 
operator uses the mouse to manually position the cursor 
at the center of the fiber. The computer 14 then 
automatically fills in the interior of the fiber and 
keeps track of the perimeter and area. Unlike Geckle 

10 and Jenkins, who employ a variation of the blob 
counting technique, an algorithm known as region 
growing is employed. 

Region growing is a rather general term. A more 
specific term for the procedure described is " flood 

15 fill" — which is usually used in the context of 
computer graphics. Region growing requires some 
starting point from where the region grows. Each pixel 
in the region expands into its 4 -neighborhood if the 
neighbor satisfies some properties. Since the axons 

20 appear to be brighter than the myelin sheaths, it can 
be assumed that the region will be bounded by dark 
pixels. The average intensity of the local 20X20 grid 
is calculated. Pixels with intensity lower than this 
average are considered to be boundary pixels. 

25 There are two ways to implement region growing. 

An easy and straight forward method is to use 
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recursion. This has been proved to be inappropriate on 
the PC because of its limited stack size. Stack 
overflows occur frequently because recursive calls fill 
up the stack rather quickly. Instead, the system 
5 according to the present invention maintains the 
locations of the outer pixels of the growing region in 
a linked list. Pixels which have been determined to be 
non-expandable, that is, those whose 4 -neighbors are 
all boundary pixels, will be deleted from the list. 
10 Expansion continues until the list becomes empty. 

Table 5.1 provides a comparison of the two 
methods. Region growing searches for the boundary 
symmetrically in ail directions while blob counting 
follows the boundary. 





i 

Region Growing 


Blob Counting . 


; Time 


0(N 7 ) 


O(N) 


! Space 


0(N 7 ) 


O(N) 


j Symmetry 


Good 


Poor 


: 

; Reliability 


iLXceuen: 


Fair 



15 This means that the computational time for region 
growing increases with the area while that of blob 
counting increases with the perimeter. This makes blob 
counting more efficient. However, when the radius of 
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che search area is small, the computational advanraoe 
cf blob counting is not that great. The average radius 
of the fibers that are counted in the present 
application of the system is only about 5 pixels. This 
5 advantage is also lessened by the fact that the 
computational time for either of these method is much 
longer than the time required by the operator to place 
the mouse at the center of each fiber. 

The reason for using region growing is the fact 
10 that the boundaries of many fibers are not nicely 
connected. The entire interior of the fiber could be 
filled before "leaking" occurs. This means that there 
is a chance to stop the operation just before the 
leaking and still have the filled area, which 
15 corresponds zo the axon area, measured. Blob counting, 
on the other hand, will not be able to handle the case 
when there are discontinuities in the boundary because 
the area cannot be measured until the whole boundary is 
traced . 

20 The AIDED mode of operation also allows the user 

to limit the search area to a grid smaller than the 
default 20X20. Limiting the grid size has two effects: 
1) it limits the search area to prevent leaking, and 2) 
it makes the calculation of the iocal average intensity 

25 more accurate. The user is also permitted to 
repetitively delete the last fiber counted if an error 
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is detected or randomly to delete any counted fiber. 
The procedure is illustrated in Figure 5.2 which shows 
a magnified region of a partially counted section. 

Problems arise when processing a fiber which 
5 touches another one. The myelin sheath area will 
sometimes grow towards the neighboring fiber. To solve 
this problem, one must restrict the growing of the 
myelin sheath area to pixels which have lower intensity 
than the average, but also only allow pixels to grow to 
10 a neighbor having a lower intensity. However, it is 
not uncommon that the traced object will still contain 
unwanted protrusions (illustrated by Fig. 5.3(a)) which 
can be eliminated by the computer performing a five 
step process of: 

15 1. Using erosion to disconnect the neighboring 

fiber. 

2. Using region growing to retrieve the area that 
is not disconnected. 

3. Using dilation to balance the effect of the 
20 first erosion. 

4 . Using dilation to smooth discontinuities in 
the fiber boundary, and 
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5. Using erosion to balance the effect of the 
second dilation. 

The fourth and the fifth steps essentially 
comprise a closing operation. For the dilation and 
5 erosion, a 3 by 3 square shape structuring element is 
used. Steps 1-5 above are respectively illustrated in 
Figs. 5. 3(b) -(f). 

5.4 AUTO MODE 

The AUTO mode of operation is in many ways similar 
10 to AIDED mode. The main and significant difference is 
that AUTO does not require operator assistance in the 
fiber identification process. In AIDED, the operator 
points out where the fibers are and this prevents the 
computer from performing an extensive search of the 
15 image in order to locate all the possible structures. 
The AUTO mode has to perform this search and most 
importantly, it has to be fast and accurate. 

The following process steps are performed by the 
computer 14 in the AUTO mode of operation: 

20 1- Identification of all the possible locations 

of the fibers. 

2. Elimination of the false identifications of 
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step 1 by thoroughly checking the neighborhood of the 
possible locations. | 

3. Measurement of the structures that have been 
determined to be fibers . 

5 5.4.1 Preprocessing 

Prior to performing the above steps, the digitized 
image should be preprocessed to produce a lowpass and 
highpass version using the image techniques described 
above. The highpass version with enhanced edges is 
10 used in the fiber detection process. The lowpass 
version with spikes removed is used in the location of 
axon and myelin sheath area. 

The image should also be corrected geometrically 
because the aspect ratio of the digitizing system is 
15 4:3. This means that a circular fiber is not reaily 
circular in the frame memory, but the horizontal 
diameter is 0.75 of the vertical diameter. The image is 
corrected by shrinking the image in the y direction: 

' .v- 3 -*vmod4) , , 1-iu mod 4) . . 

u,,/ 1 M*. U-23yJj — ;i.33yij (5.4) 

wnere ana ::r,y; are the orieinai ana corrected image rescectiveiv. 
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5.4.2 Locating Local Maxima 

The locai maxima of the image are obtained by the 
computer 14 after preprocessing is completed. Local 
maxima are pixels with higher intensity than the 
5 neighboring pixels, representing either the bright 
regions that are at the center of the fibers or the 
inter-fiber holes. it is reasoned that, since the axon 
area is brighter than the myelin sheath area, there 
must be at least one, but possibly several, local 
10 maxima located therein. 

A local maximum is defined to be a point (x,y) in 
an image i(x,y) such that the following conditions are 



met : 



> its, y _ : ; 
*l'*.y) > :(a~l T y _ i) 



The definition of local maximum insures that only one 
pixel would be identified as a local maximum even if 
there is a region of constant maximum intensity. The 
typical number cf local maxima in an image is about 
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3,000 in the present application. An example of the 
ohtained local maxima, indicated as black dots, is 
shown in Figure 5.4. 

5.4.3 Variable Thresholding 
5 Thresholding is usually successful when the image 

in question has a bimodal histogram, in which case the 
threshold value would be the intensity corresponding to 
the trough of the histogram. The histogram of the 
sample image is shown in Figure 3.1(b) which is clearly 
10 unimodal. However, thresholding can still be done 
using a variable threshold — a threshold that changes 
with location. 

It is assumed that the axon area has higher 
intensity and myelin sheath area has lower intensity 
15 than the average. In other words, the threshold value 
is the local average. Unfortunately, calculating the 
average for each point is a time consuming process. 
Instead, it is possible to calculate the average for 
221 equally spaced grid points arranged in a 17x12 
20 network. This separates the image into 132 regions of 
size A x£y. To calculate the average for any point in 
the image, the system need only to interpolate from the 
4 closest grid points whose average have already been 
determined . 
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For example, let (x Q ,y 0 ) be the point for which 
the system is trying to calculate the average; x^Xj 
are the left and right boundaries of the region 
respectively; y 1 ,y 2 are the top and bottom boundaries 
of the region, respectively; and i ., is the averaoe at 
(x.:,y k ). Then the average is: 

: = (ASz - B)6y - C&x + D 



wnere .4 = — - ~ : - l 



(5.5) 



Ay 

C = ilLZili 
= in 

<>y = v - i'i 

Note that A, B, C, and D only need to be calculated 
once for each region. 

10 The proof cf equation 5.5 is based on dividing the 

screen of the monitor into 16 x 12 regions cf size 
*xXav= 32 x 30 pixels. Altogether there are 17 x 13 
corner points, for each of which the local average 
intensity is calculated by averaging the local ^xXAY 

15 pixels. To calculate the local average intensity of 
any point on the screen, it must first be determined to 
which of the 16 x 12 regions the point belongs. Fig. 



WO 92/17848 



PCT/US92/02563 



49 



5.5 shows one of the 16 x 12 regions. The local 
average intensities at j(xi,yi), (x2^yi), (xirY2) and 
( X 2'Y2) are known, as they are one of the 17 x 13 
corner points . Both ^ and ^ are approximated by 
interpolation: 



221 - ill A* 

fi 21 - iu)6z 



A - 



,j s i: 



and therefore, 



*: - M = Ji: - «:i 

Ax 



By the same principle, but applied in the y 
direction, i is obtained by interpolation, 



• = — — — 

_ (*'?? — *!? — *n -r ?n)&rfly . (M? — <n)<fy (271 — *n )<?r 
At/ Ay Az 



I 



f*2t ~ *I3 — *?l -r Mi) f 1*13 — hi 



*2i - in 



r * - - 1 1 - * . l * I i - 

cz - ■ oy 01 -r / 

Ax Ay Av Ax 
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Thus, each of the 16 x 12 regions is characterized by 
its own A, B, C, D. 

The resulting image using the threshold calculated 
this way may still have some discontinuities in myelin 
5 sheath boundaries. These discontinuities are closed by 
a closing operation using a 2 by 2 square shaped 
structuring element. A thresholded image is shown in 
Figure 5.6. 

5.4.4 Refining Local Maxima 
10 Th © number of local maxima at this point may be 

quite large. An effective way to eliminate some of 
them is to compare the location of the local maxima 
with the thresholding image. We assume the axon areas 
to be above the threshold value and therefore the local 
15 maxima that are located at a point where the intensity 
is below the average can be eliminated. 

It should be noted that the information about the 
locai maxima is binary, that is, a point is either a 
local maximum cr not. The thresholded image is also 
20 binary, each point is either above or beiow the 
threshold. This means that the system can store the 
local maxima bit by bit in an array, -each bit 
representing one point. The bit is set to one if the 
point is a locai maxima or to zero if the point is not. 
25 a local maxima. The system stores the threshold image 
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similarly in an array, wherein a bit is 1 if the 
intensity corresponding to jthat point is above the 
. threshold and 0 otherwise. This is very space 
efficient — only 3 2 KBytes for each array. Moreover, 
5 the local maxima refinement procedure is simply a 
bitwise AND operation of the local maxima and the 
thresholded image arrays, which require less than one 
second to complete. 

5*4.5 Blob Coloring 

10 The threshold image (Figure 5.6) clearly shows 

that the image is segmented into regions. Some of 
which correspond to fiber areas and some to inter-fiber 
holes . The goal is to extract those that correspond to 
fiber areas. It is necessary, however, to locate ail 

15 of the segmented regions first. This poses a difficult 
problem because the space requirement to represent the 
segmented regions must be considered as discussed 
above . 

It is impossible to store the locations of all the 
20 points belonging to a certain region, but it is 
possible to indicate to which region a certain point 
belongs. This technique is known as blob coloring. 
The computer assigns the same color to each point that 
belongs to the same region. The color is just a value 
25 in the image frame memory. The interpretation of the 
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10 



25 



value in the image frame memory depends on the 
operation. For example, for a black and white display 
monitor, this value represents intensity or brightness. 
For blob coloring, this value represents the region to 
which the pixel belongs. The term "blob coloring" is 
used because an effective way to display the different 
regions of the image is to use a different color for 
each region. In the context of blob coloring, a color 
is just a value which is stored in the image from 
memory; a pixel having a value k belongs in region k. 
If a 12-bit image frame memory is employed, one can 
have at most 4096 different regions, one of which is 
reserved to represent the background. 

The system starts blob coloring at a local maximum 
15 because every fiber has at least one local maxima. For 
each local maximum, the region growing procedure is 
used, letting it grow up to the thresholded boundary. 
Then the system lets it grow into the myelin sheath 
area. The system then colors the region as it grows, 
using a different color for each local maximum, and 
making sure that it doesn't grow into a region that has 
already been colored. 

The resulting number of regions is always less 
than the number of local maxima. This is because the 
regions that are determined to be too big (by 
comparison to a predetermined limit) probably do not 



20 
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belong to a fiber and are discarded. Regions that have 

extremely high eccentricity (t>40) are also thrown away. 

i 

Furthermore, it is not unusual that several local 
maxima are located in the same region, in which case 
5 only one region results . 

The result of blob coloring is that many regions 
are connected to many other regions. For reasons that 
will be explained below , it is desirable to obtain 
blobs that correspond to fibers connected to inter- 
10 fiber holes, but not to other fibers. As shown in 
Figs. 5.7(a)-(b), the fiber-fiber contact area is 
smaller than fiber-hole contact area. To eliminate 

fiber-fiber contacts while preserving fiber-hole 

» 

contacts, colored pixels which have neighbors of more 
15 than two colors are deleted. This procedure is similar 
to erosion but there is no structuring element and the 
deletion is more selective. The system performs this 
"selective erosion" twice, each time possibly deleting 
one pixel from each end. 
20 If the model shown in Figs. 5.7(a)-(b) were 

correct and the fibers ' radii were 5 pixels , then 
(through geometric calculations) it could be shown that 
the fiber-hole contact length would be about 5 pixels- 
- one-sixth the circumference of a circle with radius 
25 equal to 5 pixels. This means that the fiber-hole 
contact would still exist after selective erosion, but 
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but the fiber-fiber contact would be eliminated. An 
image of a nerve section after blob coloring is shown 
in Figure 5.8. 

5.4.6 Symbolic Constraint Propagation Operation 

To eliminate blobs which represent holes , it can 
be assumed that holes usually have more neighbors than 
do fibers. This is simply because of the fact that 
holes are more irregular. In addition, since most of 
the neighbors of a fiber are holes, some of which have 
. already been eliminated in the blob coloring procedure, 
the number of neighbors next to a fiber is small. 

First, the computer sets up an adjacency list for 
each region. An adjacency list is a linked list 
containing the neighbors cf a region. For example, the 
adjacency list for region o in Fig. 5,7 contains p, q 
and r; the adjacency list for region p contains o, etc. 
Once the adjacency lists of the regions have been 
created, the blob colored, the blob colored image is no 
longer be needed. Given a region, the computer can 
look up the adjacency list to find out how many 
neighbors it has and who they are. Regions are 
eliminated based on a set of criteria. This is done 
until all the remaining regions have either zero or one 
neighbor. Regions which have zero neighbors will be 
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considered to be a fiber and will not affect the 
elimination process of the other regions in any way. 

Regions which have only one neighbor are very 
likely to be a fiber. It is desirable to delete the 
5 only neighbor of a region, because that would make that 
region a fiber. On the other hand, regions with many 
neighbors are likely to be holes . It would be even 
more desirable to delere the only neighbor of a region 
if that neighbor, in turn, has many neighbors. 
10 The elimination process uniquely utilizes 

eccentricity to eliminate blob regions that do not 
represent axons, as follows: 

1. If all of the remaining regions have 0 or 1 

* 

neighbors, then go to step 4. 

15 2. If there aren't any regions that have only one 

neighbor, delete the region with the largest number of 
neighbors. If there is a tie, delete the one with the 
highest eccentricity. Then go back to step 1. 

3. Select from the regions which are neighbors of 
20 regions which have only one neighbor, the one with the 
most number of neighbors for deletion. If there are 
several regions that fit this category, select the one 
with the highest eccentricity. Then go back to step 1. 
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4 . Two types of blobs will be obtained at the end 
of this procedure — blobs which have no neighbors, and 
pairs of blobs which are neighbors of each other. in 
the latter category are usually fiber-fiber pairs, or 
5 fiber-hole pairs, but almost never a hole-hole pairs. 
Therefore, it must be decided whether to keep both, or 
to discard one of the blobs . The decision is based on 
the eccentricities of the regions. If both of the 
regions have lower eccentricity than 14 then keep both. 
10 Otherwise, delete the blob with the higher 
eccentricity. This is the last step of the fiber 
identification process . 

It should be noted that the eccentricity of 14 was 
picked through experimentation. An eccentricity of 14 
15 appeared to give the best results. However, there is 
only a slight change in accuracy from eccentricities 
ranging from 13 to 18. 

At this point, the system identifies and counts 
the fibers by a set of local maxima which correspond to 
20 the remaining blobs. This set of local maxima is used 
as the starting points for the location of the axon and 
myelin sheath areas using region growing as described 
above with respect to the AIDED mode of operation. The 
only modification is that the average used is the one 
25 obtained from interpolation. The result is shown in 
Figure 5.9. 
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5.4.7 Time Efficiency 

Measures have beeip taken to improve the time 
required by the system to perform the AUTO mode. The 
most important one is the extensive use of the feedback 
5 look-up table of the imaging board. This increases the 
speed of several procedures like dilation , erosion and 
local maxima location. Table 5.2 shows the major steps 
of AUTO mode and the percentage of time spent on each 
procedure. The most time consuming step is blob 
10 coloring which takes ' aimost 50% of the total time 
required. Thus, a logical approach to improving the 
efficiency of system operation would be to further 
optimize the blob coloring procedure. 

5 . C Results 

15 Fig. 6.1 is a table which illustrates the results 

obtained from the system illustrated in Fig. 2.1 in the 
MANUAL, AIDED and AUTO modes of operation. The time 
required for the MANUAL mode of operation ranged from 
110 to 230 minutes depending on the skill of the 

20 operator. The time required for the AIDED mode of 
operation ranged from 16-25 minutes. The AUTO mode of 
operation, however, counted the number of fibers in the 
sample cross-section and determined axon area, myelin 
sheath thickness and eccentricity within 4.8 minutes, a 
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vast improvement over the MANUAL and AIDED modes of 
operation . 

Fig. 6.2 is a table which illustrates the analyses 
of six histological sections using the AUTO mode of 
5 operation — three from normal nerve and three from 
electrically stimulated nerves. There is an observable 
difference in axon area between normal and stimulated 
fibers; however, it is unclear where along the cochlear 
nerve the axial sections were taken, and since the size 
10 of the fibers depends on their origins in the cochlea, 
it cannot be concluded from this data that the smaller 
average axon area in the stimulated nerves is a result 
of electrical stimulation. 

Of the six basic parameters generated in the AUTO 
15 mode of operation (illustrated in Figs. 6. 3 (a) -(f) for 
the image shown in Figure 1.2), only eccentricity, 
which is unitless, does not depend on the size cf the 
objects. The fiber eccentricity of stimulated nerves 
is higher than chat of normal ones, but only on the 
20 order of 1.5% to 3.2%. This could be sufficient, 
however, to indicate damage to the stimulated nerves, 
as statistical analysis has shown that there is oniy a 
5% probability that the finding of a true average 
increase in eccentricity between 1.5% and 3.2% is 
25 incorrect. 
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The invention has been described with reference to 
certain preferred embodiments thereof. It will be 
understood, however, that variations and modifications 
can be made within the scope of the appended claims. 
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WHAT IS CLAIMED IS: 

1 . A method of determining characteristics of 
nerve fibers, said method comprising the steps of: 
generating digital image data representative of an 
image of a nerve fiber cross-section using an 
electronic camera; storing said digital image data in a 
storage unit; accessing said digital image data from 
said storage unit and locating local maxima of said 
digital image data with a control unit; performing a 
variable thresholding operation on said digital image 
data; refining said local maxima; defining a plurality 
of blob image regions using blob coloring; performing 
an elimination operation to identify which of said 
plurality of blob image regions represent nerve fibers. 

2. A method of determining characteristics of 
nerve fibers as claimed in claim 1, further comprising 
the step of calculating at least one of the axon area, 
the myelin sheath area, the fiber area, the myelin 
sheath thickness, the fiber perimeter and the fiber 
eccentricity of the blob image regions that represent 
nerve fibers. 
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3. A method of determining characteristics of 
nerve fibers as claimed in claim 1, further comprising 
the step of counting the number of blob image regions 
that represent nerve fibers. 

4 . A method of determining characteristics of 
nerve fibers as claimed in claim 1 , wherein said local 
maxima are defined as a point (x,y) in an image i(x,y) 
such that the following conditions are met: 



»'(*t*y) > 
*'(=,»; > 



t(r-l,»-l) 
»'(x,y-l) 
t(x + l,y-l) 
*(*-!»») 
*(x + l,ir) 

ii>+ l,y + 1) 



5 . A method of determining characteristics of 
nerve fibers as claimed in claim 1 , wherein said step 
of performing a variable thresholding operation 
includes the step of determining a plurality of locai 
5 average intensities corresponding to regions of said 
image . 
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6. A method of determining characteristics of' 
nerve fibers as claimed in ciaim 1, wherein said step 
of performing an elimination operation includes the 
steps of a) determining how many neighbors each of said 
blob image regions has and terminating the elimination 
operation when each of said blob image regions has only 
one or zero neighbors, b) deleting the blob image 
region with the largest number of neighbors, if there 
are no regions having only one neighbor, c) selecting 
from the blob image regions which are neighbors of blob 
image regions having only one neighbor and deleting the 
one with the most number of neighbors, and d) repeating 
steps a) , b) and c) . 

7. A method of determining characteristics of 
nerve fibers as claimed in claim 6, wherein step a) 
includes the step of identifying an eccentricity of 
said blob image regions and deleting a blob image 
region having a highest eccentricity if there are two 
blob image regions with the same number of neighbors 
and step c) includes the step of deleting the blob 
image region with the highest eccentricity if there are 
two blob image regions having the same number of 
neighbors which are neighbors of blob image regions 
having only one neighbor. 
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S . A method cf determining characteristics of 
nerve fibers as claimed in claim 7 , f urrher comprising 
the step of identifying pairs of blob image regions 
which are sole neighbors of each other and deleting the 
blob image region of each pair of blob image regions 
that has the higher eccentricity of the pair if that 
blob image region has eccentricity higher than 14 . 

9 . A method of determining characteristics of 
nerve fibers as claimed in claim 7 , f urrher comprising 
the srep of idenrifying pairs of blob image regions 
which are sole neighbors of each orher and keeping both 
blob image regions if both blob image regions have 
eccentricities lower than 14. 

10. A method of determining characteristics of 
nerve fibers as claimed in claim . 1 comprising the 
further step of obtaining an eccentricity cf said blob 
image regions and eliminating blob image regions as not 
representing axons based on eccentricity thereof. 

11. In a method cf measuring diameters of nerve 
fibers using image analysis , the improvement comprising 
the steps of: 

obtaining an eccentricity characteristic for a 
plurality of blob image regions and 



WO 92/17848 



PCI7US92/02563 



64 

eliminating blob image regions as nor representing 
nerve fibers based on the eccentricity characteristics 
thereof . 

12. The method cf claim 11 wherein said 
eliminating' step comprises the steps of determining a 
number of neighbors for each blob image region has , 
deleting £ blob image region havin-g a highest 

5 eccentricity if there are two biob image regions with 
the same number of neighbors and deleting the blob 
image region with the highest eccentricity if there are 
two blob image regions having the same number of 
neighbors which are neighbors of blob image regions 
10 having only one neighbor. 

13. The method of claim 11 comprising the step of 
using blob coloring to define said blob image regions . 
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